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A new Skyrme-like energy density suitable for studies of strongly elongated nuclei has been de- 
termined in the framework of the Hartree-Fock-Bogoliubov theory using the recently developed 
model-based, derivative-free optimization algorithm POUNDerS. A sensitivity analysis at the opti- 
mal solution has revealed the importance of states at large deformations in driving the parameter- 
ization of the functional. The good agreement with experimental data on masses and separation 
energies, achieved with the previous parameterization unedfO, is largely preserved. In addition, the 
new energy density unedfI gives a much improved description of the fission barriers in ^*''Pu and 
neighboring nuclei. 

PACS numbers: 21.60.Jz, 21.10.-k, 21.30.Fe, 21.65. Mn, 24.75. 



I. INTRODUCTION 

One of the focus areas of the UNEDF SciDAC collab- 
oration [1, 2] has been the description of the fission pro- 
cess within a self-consistent framework based on nuclear 
density functional theory (DFT). Until now, attempts 
at going beyond the macroscopic-microscopic methods 
[3] have been carried in the context of the original self- 
consistent nuclear mean-field theory [4] with Skyrme (see, 
e.g., Refs. [5-10]), Gogny (see, e.g., Refs. [11-14]), and 
relativistic (see, e.g., Refs. [5, 15, 16]) energy density 
functional (EDFs). The fundamental assumption of 
the nuclear DFT is that one can describe accurately 
a broad range of phenomena in nuclei, including ex- 
cited states and large-amplitude collective motion, by 
enriching the density dependence of the functional while 
staying at the singlc-rcfcrcncc Hartrec-Fock-Bogoliubov 
(HFB) level. In this picture, beyond-mean-field correc- 
tions are implicitly built-in. Preliminary studies aimed 
at re-examining the old problem of restoring broken sym- 
metries in this context are promising [17, 18]. 

A common challenge to both the self-consistent mean 
field and the DFT approach is the need to carefully op- 
timize the EDF parameters to the preselected pool of 
observables [4, 19-25]. In particular, special attention 
must be paid to optimize the parameters in the same 
regime where the theory will later be applied and to 
choose the fit observables accordingly. In a recent work 
[26], we showed that existing Skyrme EDFs exhibit a 
significant spread in bulk deformation properties, and 
re-emphasized [27, 28] that the resulting theoretical un- 
certainties could be greatly reduced by considering data 
corresponding to large deformation in the optimization 
process. Let us recall that the early Skyrme-type EDF 
SkM* [29] was in fact optimized for fission studies in the 



actinide region by considering the experimental informa- 
tion on the fission barrier of ^^^Pu. However, the op- 
timization was not performed directly at the deformed 
HFB level; instead, a semi-classical approach was used 
based on the Thomas- Fermi approximation together with 
shell-correction techniques. The DIS parameterization 
of the finite-range Gogny force was also fine-tuned to the 
first barrier height of ^^^Pu [30], considering a rotational 
correction to the energy of the deformed state. How- 
ever, this fine-tuning again was not done directly at the 
HFB level but by a manual readjustment of the surface 
coefficient of the EDF using a phenomenological model. 
Also, in the Bskl4 EDF of the HFB-14 mass model [31] 
by the BruxcUes-Montrcal collaboration, data on fission 
barriers were utilized to optimize the EDF parameters by 
adding phenomenological collective corrections, includ- 
ing a rotational one. One may, therefore, conclude that 
no EDF has ever been systematically optimized at the 
deformed HFB level (and without phenomenological cor- 
rections added) by explicitly considering constraints on 
states at large deformations. 

In a previous study [25], we applied modern opti- 
mization and statistical methods, together with high- 
performance computing, to carry out EDF optimization 
at the deformed HFB level, namely, the approximation 
level where the functional is later applied. The result- 
ing EDF parameterization unedfO yields good agree- 
ment with experimental masses, radii, and deformations. 
The present work represents an extension of [25] to the 
problem of fission. In particular, it builds on the results 
reported in Ref. [26], which concluded that the data on 
strongly deformed nuclear states should be considered in 
the optimization protocol to constrain the surface prop- 
erties of the functional. 

Here we propose the new EDF Skyrme parameteriza- 
tion, UNEDF 1, which is obtained by adding to the list 
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of fit observables the experimental excitation energies of 
fission isomers in the actinides. To ensure that the func- 
tional can be used in fission and fusion studies, we have 
removed the center-of-mass (c.o.m.) correction in the 
spirit of the DFT. As in the case of unedfO, a sensitiv- 
ity analysis has been performed at the solution in order to 
identify possible correlations between model parameters 
and assess the impact of the new class of fit observables 
on the resulting parameterization. 

This paper is organized as follows. In Sec. II we briefly 
review the theoretical framework, establish the notation, 
and justify the removal of the c.o.m. correction. Section 
III defines the set of fit observables, discusses numerical 
precision and implementation, and presents the new UN- 
edfI parameter set together with the results of the sen- 
sitivity analysis. To assess the impact of fission-isomcr 
data, we compare unedfI with unedfO in Sec. IV. In 
Sec. V we study the performance of unedfI with re- 
spect to global nuclear observables, spectroscopic prop- 
erties, fission, and neutron droplets. Section VI contains 
the main conclusions and lays out future work. 



details. (For brevity, we have omitted the explicit depen- 
dence of the densities on the coordinate r.) The Jt den- 
sity is the vector part of the spin-current density tensor 
J^i/. (As in our previous work [25], non- vector compo- 
nents of Jpjy were disregarded.) The coupling constants 
are real numbers, except for , which is taken to be 
density-dependent: 



(3) 



All volume coupling constants (Cf and Cf^) can be 
related to the constants characterizing the infinite nu- 
clear matter [25], and this relation was used during the 
optimization in order to define the range of parameter 
changes. 

The Coulomb contribution is treated as usual by as- 
suming a point proton charge. The exchange term was 
computed at the Slater approximation: 



cCoul 
<--Exc 



(r) 



1/3 



(4) 



II. THEORETICAL FRAMEWORK 

A. Time-Even Skyrme Energy Density Functional 

In the nuclear DFT, the total binding energy E of the 
nucleus is a functional of the one-body density p and 
pairing p matrices. In its quasi-local approximation, it 
can be written as a 3D spatial integral: 

E[p, p] = f (fr H{r) 

d?r [f^'"(r)+xo(r)+XiM 

+x(r)+f§r'(r)+fg-Hr-)] , (1) 

where ?^(r) is the energy density that is quasi-local (it 
usually depends on derivatives with respect to the local 
densities), time-even, scalar, isoscalar, and real. It is usu- 
ally broken down into the kinetic energy (f^™(r)) and 
nuclear potential (for both the particle-hole and particle- 
particle channels, Xo,i(''") and x{r), respectively) and 
Coulomb terms (fg°^'(T-) and £^°f{r)). For Skyrme 
functionals, the particle- hole energy density Xo('")+Xi('') 
splits into Xo{f)i depending only on isoscalar densities, 
and Xi('')! depending on isovector densities (and the 
isoscalar particle density through the density dependence 
of the coupling constant Cf; see below) [4, 32, 33]. Each 
term takes the generic form 



-CTptrt + cfjf 
+C^^''ptApt +CP^'ptV 



(2) 



where pt, Tt, and Jt {t = 0, 1) can all be expressed in 
terms of full density matrix pt{ra,r'a'); see Ref. [4] for 



For the pairing energy density x(^)? "we use the mixed 
pairing description of [34] with 



X(r) = E 



•^0 



1 - 



1 Po{r) 



Pc 



P'ir), 



(5) 



where p is the local pairing density. The value pc=0.16 
fm~^ is used throughout this paper. We allow for differ- 
ent pairing strengths for protons (Vg') and neutrons {V^) 
[35]. A cut-off of -Ecut = 60MeV was used to truncate 
the quasi-particle space [36]. To prevent the collapse of 
pairing, we used the Lipkin-Nogami procedure according 
to [37]. 



B. Treatment of the Center of Mass 

The success of the self-consistent mean-field theory is, 
to a great extent, due to the concept of symmetry break- 
ing. A classic example is the breaking of the translational 
invariance by the mean field that is localized in space. 
The associated c.o.m. correction to the binding energy 
[4, 38, 39], -{Pl^ )/{2mA), is usually added to the DFT 
binding energy in (1). This correction contributes typ- 
ically a few MeV to the total energy. Moreover, it was 
shown that adopting approximations to this correction 
during the optimization of the functional could lead to 
significantly different surface properties [4, 17, 40]. 

Since the c.o.m. correction is not additive in particle 
number, it causes serious conceptual problems when deal- 
ing with fission or heavy-ion fusion, that is, when one 
considers the split of the nucleus into several fragments, 
or formation of the compound nucleus through a merger 
of two lighter ions. In fission studies, it was shown that 
the contribution of the c.o.m. correction between the two 
prcfragmcnts could amount to several MeV near scission 
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[41-43]. Moreover, properly computing this relative con- 
tribution is difficult, as it reflects the degree of entan- 
glement between prefragments [43, 44]. Time-dependent 
Hartree-Fock calculations of low-energy heavy-ion reac- 
tions are even more problematic, as there is currently no 
solution to the discontinuity of the c.o.m. correction be- 
tween the target -f-projectilc system and the compound 
nucleus [45, 46]: such calculations usually neglect the 
c.o.m. term altogether, even though EDFs employed in 
such calculations have been usually fitted with the c.o.m. 
correction included. There are, however, some excep- 
tions, see, e.g., Ref. [47]. Note that the same problem oc- 
curs with the so-called rotational correction arising from 
the breaking of rotational invariance by deformed mean 
fields [41, 48]. 

Another undesired property of the c.o.m. correction is 
that it slightly breaks the variational nature of HFB when 
adding or subtracting a particle [49], i.e., it violates the 
Koopmans theorem. The resulting s.p. energy shifts are 
quite significant and they are of the order of the mass 
polarization effect related to the fact that when adding 
or subtracting a particle to a closed spherical core, the 
resulting nucleus becomes deformed. If a spherical sym- 
metry is imposed, the mass polarization effect is a self- 
consistent rearrangement of all nucleons, when an odd 
particle is introduced to the system. This corresponding 
energy shift is E^i^/A^, i.e., about 0.4 MeV in "^OCa [49]. 

As discussed in Sec. I, the EDF is supposed to cap- 
ture all the physics of interest at the HFB level. In 
other words, the functional is to be built from the full 
single-reference density matrix. While the HFB vacuum 
breaks the translational symmetry, the associated correc- 
tion term should be absorbed in the density dependence. 
In particular, it should be possible to express the c.o.m. 
term not as an explicit function of A, as currently being 
done, but through a density functional. Until a simple 
prescription is proposed, however, it is more consistent to 
simply drop all corrections that originate from a Hamil- 
tonian view of the problem, including the c.o.m. term. 
And, since our focus in on fission, this is precisely what 
we have done in this work. 



TABLE I: Experimental excitation energies of fission isomers 
[50] (in MeV) considered in the unedfI dataset. 
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N 


E 


92 


144 


2.750 


92 


146 


2.557 


94 


146 


2.800 


96 


146 


1.900 



and four excitation energies of FIs. For the FIs, we used 
the weight Wi = 0.5 MeV in the objective function 

j \ * / 

The weights for binding energies, proton rms radii, 
and odd-even mass (OEM) staggering are the same as in 
Ref. [25]. 

Two assumptions made in [25] were also adopted here: 
(i) since the isovcctor effective mass cannot be reliably 
constrained by the current data, it was set to 1/A/* = 
1.249 as in unedfO and the SLy4 parameterization [51], 
which was the initial starting point in our optimization; 
and (ii) since tensor terms are mostly sensitive to the 
single-particle (s.p.) shell structure, which is not directly 
constrained by the unedfI dataset. the tensor coupling 
constants Cg and C'l were set to zero. In summary, 
compared with unedfO [25], the optimization of unedfI 
is characterized by the following: 

• The same 12 EDF parameters to be optimized, 
namely, p,, E^^/A, K^^, <m, L^m^ ^j. ^ ^^p ^ 

Cf'. Vl\ V^. CS^', and Cf^^; 

• 7 additional data points: 3 new masses and 4 FI 
energies with the weights w = 0.5 MeV; 

• Neglect of the c.o.m. correction term. 



III. OPTIMIZATION AND SENSITIVITY 
ANALYSIS 

A. Experimental Dataset 

Since the focus of this study is the construction of 
an EDF optimized for fission, our experimental dataset 
has been expanded by including the excitation energies 
(bandheads) of four fission isomers (FIs) listed in Ta- 
ble I. The ground-state (g.s.) binding energies of ^^^U, 
^^°Pu, and ^"^^Cm, not included in the previous unedfO 
dataset [25], were added for consistency (the g.s. energy 
of 236 u Tj^g^g already there). Consequently, compared with 
unedfO, the unedfI dataset contains seven new data 
points: three additional g.s. masses of deformed nuclei. 



B. Numerical Precision and Implementation 

All HFB calculations were run with the code hfbtho 
[52]. The code expands the HFB solutions on the axially 
symmetric, deformed harmonic oscillator (HO) basis. In 
the optimization of unedfO, we used a spherical basis 
with 20 HO shells, which was found to give a good com- 
promise between the numerical precision and computa- 
tional performance. The current optimization includes 
states with much larger deformation than in the ground 
state, and the dependence of the energies with respect to 
the set of basis states is more significant. 

In the unedfI optimization, all quantities but the four 
fission isomers were computed with the spherical HO ba- 
sis of Nsh = 20 shells, which includes N = 1771 basis 
states. For the fission isomers, we adopted a stretched 
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HO basis with deformation /3 = 0.4. The basis con- 
tains up to Nsii = 50 oscillator shells with an upper 
limit of iV = 1771 basis states with lowest HO s.p. en- 
ergies. The oscillator frequency uJq = oj^oj\\ was set at 
HuJo = 1.2 X 41/^1/3 jYjgV. As seen in Fig. 1, at this se- 
lection of the HO basis, the dependence of FI energies 
on the basis deformation remains fairly constant around 
/3 = 0.4. Moreover, the range of variations is significantly 
less than the corresponding weight, Wi — 0.5 McV. 




0.0 0.2 0.4 0.6 



Basis Deformation (32 

FIG. 1: (Color online) Excitation energies of fission isomers 
considered in the unedfI optimization as functions of the HO 
basis deformation. 

Optimization calculations were performed on Ar- 
gonne National Laboratory's Fusion cluster, managed 
by Argonnc's Laboratory Computing Resource Center 
(LCRC). Fusion consists of 320 computing nodes, each 
with dual quad-core Pentium Xeon processors. By us- 
ing Intel's Math Kernel Library and the Intel Fortran 
compiler (ifort), we were able to run hfbtho in almost 
half the time when compared with the prebuilt reference 
BLAS library implementation and GNU's gfortran com- 
piler. We were also able to dramatically reduce the wall- 
clock time of an hfbtho computation by using OpcnMP 
at the node level to parallelize key computational bottle- 
necks. These bottlenecks involved iteratively computing 
the eigenvalues and eigenvectors of the (Jl, tt) blocks of 
the HFB matrix, as well as the density calculations re- 
flecting the same block pattern. OpenMP allowed us to 
dynamically assign processors to blocks of data for paral- 
lel processing, which further reduced the wall-clock time 
by a factor of 6 when running on an eight-core node. 

The parameter estimation computations presented in 
this paper ran 218 total simulations of hfbtho for each 
nucleus in the dataset, using 80 compute nodes (640 
cores) for 5.67 hours. As highlighted in [25], using the 
POUNDerS algorithm (Practical Optimization Using No 
Derivatives (for Squares)) on the type of fitting prob- 
lem considered here requires more than 10 times fewer 



HFBTHO runs over a more traditional, derivative-free 
Nelder-Mead optimization method [53]. Hence, with- 
out the algorithmic and computational advancements de- 
tailed above, a similar optimization could have previously 
consumed a month of computations using 80 cores of the 
Fusion cluster. 

We emphasize that, strictly speaking, both the UN- 
edfO and the unedfI parameterizations obtained in this 
work should always be used in their original environment. 
In particular, the pairing EDF should be that of Eq. (5) 
used with the original pairing space cut off; pairing calcu- 
lations must be complemented by the Lipkin-Nogami pre- 
scription; and the proton and neutron pairing strengths 
must not vary from the values determined by our opti- 
mization. In short, contrary to usual practice, there is 
no flexibility in the treatment of the pairing channel. 



C. Result of the Optimization: UNEDFI 
Parameter Set 

The starting point for our POUNDerS optimization was 
the previously obtained unedfO parameterization. After 
177 simulations, the algorithm reached the new optimal 
result. The resulting parameter set is listed in Table II. 
The first six parameters were restricted to lie within finite 
bounds, also listed in Table II, that were not allowed to 
be violated during the optimization procedure. As can be 
seen, parameters E^^/A and are on the boundary 

value. In the case of UNEDfO, we recall that and 
1/M* also ended up at their respective boundaries. The 
saturation density pc is given with more digits than the 
other parameters. Such extra precision is needed when 
computing volume coupling constants [25]. 



TABLE II: Optimized parameter set unedfI. Listed are 
bounds used in the optimization, final optimized parameter 
values, standard deviations, and 95% confidence intervals. 



X 


Bounds 




(7 


95% CI 


Pc 


[0.15,0.17] 


0.15871 


0.00042 


[ 0.158, 0.159] 




[-16.2,-15.8] 


-15.800 






^-NM 


[220, 260] 


220.000 






NM 


[28, 36] 


28.987 


0.604 


[ 28.152, 29.822] 


r NM 


[40, 100] 


40.005 


13.136 


[ 21.841, 58.168] 


1/m: 


[0.9, 1.5] 


0.992 


0.123 


[ 0.823, 1.162] 




[— cxj, +co] 


-45.135 


5.361 


[ -52.548, -37.722] 


QP^P 


[— cxj, +co] 


-145.382 


52.169 


[-217.515, -73.250] 




[— cxj, +co] 


-186.065 


18.516 


[-211.666,-160.464] 


Vi 


[— oo, +oo] 


-206.580 


13.049 


[-224.622,-188.538] 


r^p^J 

'-'0 


[— oo, -f oo] 


-74.026 


5.048 


[ -81.006, -67.046] 


(jpVJ 


[— oo, +oo] 


-35.658 


23.147 


[ -67.663, -3.654] 



We first note that the same minimum was obtained 
by starting either from the unedfO solution or from the 
UNEDFlex parameterization discussed below: this gives 
us confidence that the parameter set listed in Table II 
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is sufficiently robust. We can tlien observe tfiat most of 
the parameter values of unedfI are fairly close to those 
of unedfO [25]. There are, nevertheless, a couple of no- 
table exceptions. First, the magnitude of Cf^'' is now 
much larger. This is potentially dangerous, as it might 
trigger scalar-isovector instabilities in the functional that 
could appear in neutron- rich nuclei [54, 55]. (Our mass- 
table calculations with unedfI do not show indications 
of instability in even-even nuclei.) Second, has 
drifted considerably from its initial value, even changing 
sign. These two coupling constants control the isovector 
surface properties of the nucleus; hence, only proper con- 
straints on the shell structure like, for example, spin-orbit 
splitting in neutron-rich nuclei will allow these terms to 
be pinned down. For the moment, both coupling con- 
stants are relatively unconstrained, as evidenced also by 
their relatively large cr-value shown in Table II. 

D. Sensitivity Analysis 

1. Correlation Matrix of unedfI 

We have performed a sensitivity analysis at the solu- 
tion of the optimization. All residual derivatives were 
estimated by using the optimal finite-difference proce- 
dure detailed in Rcf. [56]. Since some of the parameters 
ran at their bounds during the optimization, the sensitiv- 
ity analysis was carried out in a subspace that does not 
contain these parameters. The same strategy was also 
used in the previous sensitivity analysis of the unedfO 
parameterization; we refer to [25] for a detailed discus- 
sion of the available options in the case of constrained 
optimization. In Table II we list the standard deviation, 
a, and 95% confidence interval (CI) for each parameter 
at the solution. As discussed in Sec. IIIC, the standard 
deviations of most of the parameters is relatively small. 

The correlation matrix for the unedfI parameter set 
is presented in Table III. It was calculated as in [25] 
and, similarly, corresponds to the 10-dimensional sub- 
space of the parameters that are not at their boundary 
value. Generally, most of the parameters are only slightly 
correlated to each other, with a few notable exceptions 
(correlations below 0.8 are not very significant from a sta- 
tistical viewpoint). The strong correlation between 1 /M* 
and both Vj^" and Vq had already been noticed in the UN- 
EDFnb parameter set of Ref. [25] and reflects the inter- 
play between the level density at the Fermi-surface and 
the size of pairing correlations. Similarly, both pairing 
strengths arc strongly correlated with Cq'^'', which can 
also be related to surface properties of the functional. In- 
terestingly, both pairing strengths are now strongly cor- 
related with one another, which was not the case with 
unedfO. The same correlation matrix of Table III is 
shown graphically in Fig. 2. 

Next, we study the overall impact of each data type in 
our function on the obtained parameter set. As in [25] 
we calculate the partial sums of the sensitivity matrix for 




FIG. 2: (Color online) Absolute values of the correlation ma- 
trix of Table III presented in a color-coded format. 

each data type. Let us recall that the sensitivity matrix 
S is defined as 

5(x)= [,/(x)jT(x)]-\/(x), (7) 

where J(x) is the Jacobian matrix. The results are illus- 
trated in Fig. 3, where we have summed absolute values 
of each data type for each parameter. The total strengths 
for each parameter were then normalized to 100%. Note 
that the fission isomer excitation energies represent less 
than 4% of the total number of data points but account 
for typically 30% of the variation of the parameter set. In 
the case of the symmetry energy coefficient, this percent- 
age is even 75% (see Sec. IV for more discussion). Com- 
pared with unedfO, we find that the overall dependence 
on the proton radii has significantly decreased, except for 
p and Ci^'' , and that the dependence on the OES has 
actually increased. This kind of analysis, however, docs 
not address the importance of a particular data point to 
the obtained optimal solution. 

A complementary way to study the impact of an indi- 
vidual datum on the obtained parameter set is therefore 
presented in Fig. 4. Here, we have plotted the amount 
of variation 




for the optimal parameter set when data points di are 
changed by an amount of O.lwi one by one. As can be 
seen, the variations are small overall, assuring us that the 
datasct was chosen correctly. The masses of the double 
magic nuclei ^°^Pb and ^'^Ni seem to have the biggest rel- 
ative impact on the optimal parameter set. One can also 
see that the sensitivity of the parameters on the new FI 
data is larger than the average datum point. By contrast, 
the dependence of the parameterization on the masses of 
deformed actinidc and rare earth nuclei is weaker. 
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TABLE III: Correlation matrix for unedfI parameter set. Correlations greater than 0.8 (in absolute value) are in boldface. 



Pc 


1.00 




















NM 


-0.35 


1.00 


















T NM 

-^sym 


-0.14 


0.71 


1.00 
















1/71 /T* 


O.o2 


U.2o 


U.OD 


1.00 
















-0.25 


-0.25 


-0.35 


-0.99 


1.00 














-0.06 


-0.15 


-0.77 


-0.22 


0.19 


1.00 












-0.32 


-0.22 


-0.36 


-0.99 


0.98 


0.22 


1.00 










-0.33 


-0.18 


-0.29 


-0.97 


0.97 


0.15 


0.96 


1.00 








-0.14 


-0.20 


-0.32 


-0.86 


0.91 


0.22 


0.85 


0.84 


1.00 






0.05 
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-0.13 


-0.10 


0.07 
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0.10 


0.07 
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FIG. 3: (Color online) Sensitivity of unedfI to different types 
of data entering the function. 



2. Discussion of the Coulomb Exchange Term 



FIG. 4: (Color online) Overall change in x for the unedfI 
when the datum di is changed by an amount of 0.1 one by 
one. The four rightmost data points marked FI correspond to 
excitation energies of fission isomers. The results for unedfO 
and UNEDFlex of Sec. HID 2 are also shown. 



It has been argued in [57] that removing the Coulomb 
exchange term from the functional could improve the 
overall fit on nuclear binding energies. This procedure 
had been motivated by the earlier works Refs. [58, 59], in 
which similar ameliorations, albeit on a smaller dataset, 
were observed. The origin of these ad hoc manipulations 
was the observation that many-body effects induced by 
the long-range Coulomb force among protons manifest 
themselves in the form of a (positive) correlation energy, 
which, to some extent, can cancel out the (negative) ex- 
change term [60, 61]. Since such an exchange-correlation 
effect is absent from the standard Skyrme functional, one 
could feel justified to simulate it by effectively screen- 
ing the Coulomb exchange term with an empirical factor 
< Qfex < 1- The special cases aox = and aex = 1 give, 
respectively, the case without and with full Coulomb ex- 
change. 

The result of the optimization of the functional with 



this additional parameter ctox is given in Table IV. Our 
objective function is slightly decreased from 51.058 to 
49.341 when this term is present. Overall, both param- 
eterizations, with and without the Coulomb exchange 
screening term, are very similar. However, one can see 
that the 95% CI is relatively large for the screening 
parameter, the value of which is also close to 1 (full 
Coulomb exchange). We recall that this confidence in- 
terval is extracted from the correlation matrix computed 
in the 10-D space of "inactive" parameters, namely, the 
space of the 10 parameters that are not at their bound 
and thus actively constrained. If one computes the Jaco- 
bian matrix in the original 13-D space of all parameters 
with a tangent plane approximation to account for the 
three active parameters E^^, K^^, and Lf^, we find 
that the 95% CI for the screening parameter becomes 
[-1.663,3.290]. This implies that ctox is basically not con- 
strained with the current dataset. 
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TABLE IV: Optimized parameter set UNEDFlex. Listed are 
bounds used in the optimization, final optimized parameters, 
standard deviations, and 95% confidence intervals. 



X 


Bounds 


X(fln.) 


a 


95% CI 


Pc 


[0.15, 0.17] 


0.15837 


0.00049 


[ 0.158, 0.159] 




[-16.2, -15.8] 


-15.800 


- 


- 


^NM 


[220, 260] 


220.000 


- 


- 


NM 


[28, 36] 


28.384 


0.711 


[ 27.417, 29.351] 


T NM 

-^sym 


[40, 100] 


40.000 






l/M* 


[0.9, 1.5] 


1.002 


0.123 


[ 0.835, 1.169] 




[— oo, +oo] 


-44.602 


5.349 


[ -51.872, -37.331] 




[— oo, +oo] 


-180.956 


47.890 


[-246.050,-115.863] 




[— oo, -f oo] 


-187.469 


18.525 


[-212.649,-162.288] 




[— oo, +oo] 


-207.209 


13.106 


[-225.024,-189.395] 




[— oo, +oo] 


-74.339 


5.187 


[ -81.389, -67.289] 




[— oo, -l-oo] 


-38.837 


23.435 


[ -70.690, -6.984] 




[0, 1] 


0.813 


0.154 


[ 0.604, 1.023] 



The dependence of every parameter on the four types 
of data included in the dataset (masses, charge radii, 
OES, and FI data) reveal an interesting consequence of 
the screening of the Coulomb exchange term. Figure 5 
shows the analogue of Fig. 3 when the screening pa- 
rameter aex is included. Note the striking difference in 
the bar plot for the symmetry energy parameter a^-^. 
fluctuations in this parameter under a variation of the 
excitation energy of the fission isomers are reduced to 
less than 10%, compared with nearly 75% when the full 
Coulomb exchange term is computed. 




FIG. 5: (Color online) Similar as in Fig. 3 but for the UN- 
EDFlex parameter set. 

This behavior can be qualitatively understood by re- 
calling a few simple facts about the bulk nuclear energy 
and deformation. A variation of the excitation energy of 
very deformed states such as fission isomers essentially 
affects the bulk surface properties of the functional — 



especially if the coupling constants driving shell effects 
are somewhat constrained by the dataset. In the lan- 
guage of the leptodermous expansion of Sec. IV C, this 
implies that both the surface and surface-symmetry en- 
ergy coefficients (which depend in a nontrivial way on the 
coupling constants of the functional) should be impacted. 
On the other hand, wc may assume that the isospin de- 
pendence of the binding energy (i.e., the total symmetry 
energy) is relatively well constrained by the several long 
isotopic sequences present in our dataset. We therefore 
see that the requirement of having the full symmetry en- 
ergy constrained together with a relatively large varia- 
tion of the surface terms should lead to a relatively large 
variation of the volume symmetry a^^, which is indeed 
observed in Fig. 4. 

One can now understand the difference of behavior of 
^sym under a change of data when the Coulomb screen- 
ing term is present: according to [60, 61], the many-body 
Coulomb correlation energy that is simulated by a^x < 1 
essentially represents a proton surface effect. Changes 
in bulk surface properties triggered by variations in the 
excitation energy of fission isomers can be entirely ab- 
sorbed by a readjustment of aex, especially since the lat- 
ter is poorly constrained by the other data, rather than 
by af^. Lastly, we note that the Coulomb exchange 
term, which is approximated by the usual local Slater 
expression, may get worse at large deformations [62]. 

In summary, considering that (i) ckox is poorly con- 
strained by the data, yet may affect significantly other 
parameters like and (ii) ctox does not significantly 
improve the quality of the fit, we decided to retain the 
full Coulomb exchange term in the present unedfI pa- 
rameterization. 



IV. CHARACTERIZATION OF UNEDFI 
PARAMETERIZATION 

In this section, we discuss general properties of the 
UNEDFl parameterization and compare it with unedfO. 

A. Energy Density in {t, x) Parameterization 

For practical applications, it is useful to express the 
coupling constants of unedfO and unedfI in the tra- 
ditional (i, a;)-parameterization of the standard Skyrme 
force, see Appendix A of [4]. The results are given in 
Table V. As can be seen, in the (t, a;)-parameterization 
the two functionals are quite different. This is to be 
expected as the relation between the C and (i, x) param- 
eterizations is partially nonlinear [63]. 

B. Energy Density Parameters in Natural Units 

The EDF parameters can also be expressed in terms of 
natural units [63] . In Table VI we list the parameter set 
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TABLE V: Parameters {t,x) of unedfO and unedfI. 



Par. 


unedfO 


UNEDFl 


Units 


to 


-1883.68781034 


-2078.32802326 


MeV-fm-' 


ti 


277.50021224 


239.40081204 


MeV-fm^ 


t2 


608.43090559 


1575.11954190 




t-i 


13901.94834463 


14263.64624708 


MeV-fm3+^^ 


Xq 


u.uuy * 4:0 * 






x-i 


-1.77784395 


-5.07723238 




X2 


-1.67699035 


-1.36650561 




X3 


-0.38079041 


-0.16249117 




64 


125.16100000 


38.36807206 


MeV-fm-"^ 


61 


-91.2604000 


71.31652223 


MeV-fm-"^ 


7 


0.32195599 


0.27001801 





of UNEDFl in standard units and in natural units. Here 
we have used the same value for the scale A = 687 MeV, 
characterizing the breakdown of the chiral effective the- 
ory, which has been found in Ref. [63] . From the numbers 
in Table VI one can see that most of the UNEDFl param- 
eters are natural, with only two minor exceptions. First, 
because the effective mass M* in our optimum is close to 
unity, the Cg^ is abnormally small. Second, Cf'^'' seems 
be on the borderline of being unnaturally large. As Table 
II indicates, however, the standard deviation for this pa- 
rameter is rather large. It has to be noted, however, that 
there is nothing unusual about the magnitude of Cf'^''. 
Indeed; some examples of EDF parameterizations with 



similar or larger values of 
[63]. 



can be found in Fig. 2 of 



TABLE VI: Coupling constants of unedfO and unedfI in 
normal units and in natural units. Value A = 687 MeV was 
used. 



unedfO 



UNEDFl 



Coupling 


Normal 


Natural 


Normal 


Natural 


Constant 


Units 


Units 


Units 


Units 


r<pp 


-706.38 


-0.795 


-779.37 


-0.878 


r'PP 


240.26 


0.271 


288.01 


0.324 


r<pp 


868.87 


0.901 


891.48 


0.937 


r<pp 


-69.77 


-0.072 


-201.37 


-0.212 




-12.92 


-0.176 


-0.99 


-0.014 


cr 


-45.08 


-0.616 


-33.52 


-0.458 


QP^P 


-55.26 


-0.755 


-45.14 


-0.616 


QP^P 


-55.62 


-0.759 


-145.38 


-1.985 




-79.53 


-1.086 


-74.03 


-1.011 


^PVJ 


45.63 


0.623 


-35.66 


-0.487 


7 


0.3220 




0.2700 





C. Leptodermous Expansion 

To extract global properties of the energy functional 
and relate them to the familiar constants of the liquid 



drop model (LDM), one needs to carry out the leptoder- 
mous expansion. The general strategy behind the expan- 
sion of nuclear EDF has been discussed in Ref. [64] , where 
one can find the relevant literature on this topic. The 
starting point is the LDM binding energy per nucleon 
expanded in the inverse radius (oc A~^/^) and neutron 
excess I ^ {N - Z)/A: 



£{A,I) = flv 



)1 + as 



(9) 



"sym-' • 



For any functional, our approach consists of combin- 
ing nuclear matter (NM) calculations with Hartree-Fock 
(HF) calculations for a large set of spherical nuclei to ex- 
tract by linear regression the various parameters of the 
expansion (9) according to the following procedure. 



TABLE VII: Liquid drop coefficients of unedfO and unedfI 
(all in MeV). 





flvol 




(2) 
^sym 






unedfO 


-16.056 


30.543 


4.418 


18.7 7.1 


-44 


UNEDFl 


-15.800 


28.987 


3.637 


16.7 8.8 


-29 



First, the bulk parameters Ovoi and agym are directly 
obtained from NM calculations. Second, the smooth en- 
ergy per nucleon £{A,I) is extracted from the spherical 
HF calculations of (A, I) nuclei by removing the shell 
correction [64]. The isoscalar coefficients of the expan- 
sion (9) can then be deduced from the smooth energy by 
plotting 

Ovol] A^/^ > a^urf + Ocurv^"^/^ 



[£{A,Q) 



(10) 



as a function of A ^Z'^. The asurf coefficient is obtained as 
the extrapolation of the curve to A-^/^ — > 0. The cur- 
vature coefficient Ocurv is then estimated from the slope 
of the line. 

The determination of isovector coefficients starts with 
the second-order symmetry coefficient Osym. It is easily 
estimated by systematic calculations in asymmetric NM. 
Defining 

a(^S(oo,/) = [£^^i(oo,/)-£:N^(oo,0)]//2 

(2) j2 y^^i 

^ f^sym I ^sym-' ; 

one can extract the second-order symmetry coefficient 
from the slope of asym(oo,/) vs. I^. Extracting the 
surface-symmetry coefficient is more involved. We first 
introduce the effective symmetry coefficient for a finite 
nucleus as 



Jeff) 



[A, I) = [E{AJ)-E{A,0)\/P 



(12) 



In nuclear matter (A — >■ -l-oo), the effective symmetry 
coefficient reduces to (11). To avoid multidimensional 
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UNEDFO 



UNEDF1 



> 




CD 


-42 








-44 


E 




>^ 
(fi 
(fi 


-46 


CO 






194 

0) 



t 

0) 

CO 



19.0 




0.04 0.08 



-28 
-30 
-32 

17.4 
17.2 
16.8 





A-1/3 




0.04 0.08 




EDF was optimized primarily to nuclear masses. Indeed, 
the main difference between unedfO and unedfI lies 
in surface properties. Relatively low values of Ogurf and 
flssym of unedfI reflect the new constraints on the FI 
data and the neglect of the c.o.m. term. Again, com- 
paring the LDM values of unedfI with those in Ta- 
ble I of Ref. [64], we note that the LDM parameters 
of unedfI are closest to those of the BSk6 EDF [651 



17.3 MeV and a.. 



-33McV) and the LSD 



LDM [66] {asnrf = 17.0 MeV and assym -38.9 MeV). 
In Sec. VC, we shall see that the reduced surface en- 
ergy of unedfI with respect to unedfO has profound 
consequences for the description of fission barriers in the 
actinides. To see this reduction more clearly, we inspect 
the effective surface coefhcient 



FIG. 6: (Color online) Surface-symmetry coefficient (15) (up- 
per panels) and surface coefficient (lower panels) versus A^^^^ 
for unedfO (left) and unedfI (right). 



regression analysis, we introduce the reduced symmetry 
coefficient by subtracting the /^-dependent part of the 
NM limit to a.. 



,(cff) 



(cff). 
sym ■ 



(rcd)^ 



ai-'iAJ)^[£{A,I)-£{A,0)] /I' 



sym 



(13) 



At the perfect LDM limit, the quantity aiy'm'iA,!) 
should not depend on the neutron excess. At small 
isospins, however, numerical uncertainties in the shell- 
correction procedure are amplified by the denomina- 
tor. In practice, it is more efficient to build an /-averaged 
reduced asymmetry coefficient. 



„(rod) 



1 



5 — a 



d/a(-^)(A,/), 



(14) 



where we choose a = 0.1 and b — 0.2 [64]. The surface- 
symmetry energy is then obtained from 



(red) 
^sym 



^1/3 



(15) 



Figure 6 illustrates the numerical accuracy of the 
method of evaluation for the surface and surface- 
symmetry coefficients of the LDM. The dashed blue lines 
indicate the fitting lines from which the final values of 
Osurf, Ocurv, and assym SLTB dcduccd. The case of unedfI 
seems to be clear. The trend of the surface energy for 
unedfO is less clean. The two groups of nuclei, huge 
and large, seem to follow slightly different slopes, and 
the fit represents a compromise. The resulting surface 
and curvature energy have to be taken with care. 

The LDM parameters of unedfO and unedfI are 
given in Table VII. As seen in Table I of Ref. [64] and 
Fig. 1 of Ref. [26], symmetry coefficients of phenomeno- 
logical LDM mass models cluster aromid cisym ~ 

30 MeV 

and a. 



-45 MeV, and the unedfO values are right 
in the middle. This result is not surprising, as this 



ssym 



^surf 



For 240pu^ the value of a'f^J is 15.33 McV for unedfI, 
16.63 MeV for unedfO, 15.87 MeV for SLy4, 15.75 MeV 
for BSk6, 15.15 MeV for SkM*, and 15.17MeV for LSD. 



flss 



■ssym 



I' 



(16) 



V. PERFORMANCE OF UNEDFI 
A. Global Mass Table 

One of the key elements required from the universal 
EDF is the ability to predict global nuclear properties, 
such as masses, radii, and deformations, across the nu- 
clear chart, from drip line to drip line. We have cal- 
culated the g.s. mass table with unedfI for even-even 
nuclei with N,Z > 8. Table VIII contains the rms de- 
viations from experiment for binding energies, separa- 
tion energies, averaged three-point odd-even mass differ- 
ences, and proton radii. Since the set of fit observables 
constraining unedfI is biased toward heavy nuclei, we 
also show rms deviations for light {A < 80) and heavy 
{A > 80) subsets. 

Figure 7 displays the binding energy residuals (i.e., 
deviations from experiment). From this figure and Ta- 
ble VIII, we can see a couple of trends. First, the en- 
ergy residuals with UNEDfI are larger than those for un- 
edfO. This result is not surprising, as the new data on FI 
and the removal of the center-of-mass correction strongly 
disfavors the lightest nuclei during the optimization pro- 
cess. Second, the characteristic arc-like behavior between 
the magic numbers is pronounced, although this trend is 
much weaker than, for example, for the SLy4 functional 
(see Fig. 7 of Ref. [25]). 

In Fig. 8 we display the residuals of two-neutron and 
two-proton separation energies. Again, the emphasis of 
unedfI on heavy nuclei is clearly seen, and the corre- 
sponding rms deviations in Table VIII quantify this fea- 
ture. Notice that two-proton separation energies are sys- 
tematically overestimated. The same trend is observed 
for the unedfO functional. We can speculate about 
sources for this effect: (i) Following the arguments of 



TABLE VIII: Root-mean-square deviations from the exper- 
imental values for unedfO and unedfI for different observ- 
ables calculated in even-even systems: binding energy E, two- 
neutron separation energy 5*2^, two-proton separation energy 
S2p, three-point odd-even mass difference A^f' (all in MeV), 
and rms proton radii _Rp (in fm). Columns are observable, 
RMS deviation for unedfO and unedfI, and the number of 
data points. 



Observable 


unedfO 


unedfI 


# 


E 


1.428 


1.912 


555 










e {a> 80) 


1.200 


1.705 


442 




0.758 


0.752 


500 


S2n {A < 80) 


1.447 


1.161 


99 


San {A > 80) 


0.446 


0.609 


401 




0.862 


0.791 


477 


Sap {A < 80) 


1.496 


1.264 


96 


Sap {A > 80) 


0.605 


0.618 


381 


^(3) 


0.355 


0.358 


442 


A^f' {A < 80) 


0.401 


0.388 


89 


Ajf' {A > 80) 


0.342 


0.350 


353 


A(^' 


0.258 


0.261 


395 


A^^' {A < 80) 


0.346 


0.304 


83 


A^^' (A > 80) 


0.229 


0.248 


312 


Rp 


0.017 


0.017 


49 


Rp {A < 80) 


0.022 


0.019 


16 


Rp {A > 80) 


0.013 


0.015 


33 



[60] , one may argue that the standard Skyrme function- 
als, such as unedfO and unedfI, lack the capability to 
describe many-body Coulomb effects; (ii) The explicit 
contribution of the Coulomb field to the pairing channel 
[67, 68] is not taken into account. It is expected that 
separate pairing strengths for neutrons and protons, as 
in unedfI, will partly account for this missing contri- 
bution [35]. However, the existence of nontrivial correla- 
tions between pairing strengths and other parameters of 
the functional (see Table III) may have consequences for 
observables such as two-proton separation energies. 

To compare unedfO and unedfI quantitatively, we 
can assess their performance on various observables listed 
in Table VIII. It is expected that since new constraints on 
fission isomers have been added when optimizing unedfI 
while keeping the same number of parameters optimized, 
the rms deviations for masses and separation energies 
must increase. Indeed, the rms deviation for the masses 
is slightly worse for UNEDfI, for both light and heavy 
nuclei. Interestingly, the quality of values remains 
roughly the same in both cases, as is true also for odd- 
even mass differences and proton radii. 
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FIG. 7: (Color online) Binding energy residuals between un- 
edfI results and experiment for 555 even-even nuclei. Iso- 
topic chains of a given element are connected by lines. 



B. Spherical Shell Structure 

The nuclear shell structure has a substantial impact on 
many nuclear properties. Notably, the single-particle lev- 
els close to the Fermi surface affect many nuclear prop- 
erties such as the strength of pairing correlations and 
dcformability. Compared with our previous work [25], 
the s.p. energies that we report here have been obtained 
from proper blocking calculations at the HFB-I-LN level 
[69] , instead of being the eigenvalues of the HF Hamil- 
tonian. This choice is motivated by the need to stay 
within a logically consistent framework: both UNEDfO 
and unedfI have been optimized at the HFB-I-LN level, 
and hence should be employed exclusively in this context. 
Moreover, in the nuclear mean-field theory with effective 
interactions, HF eigenvalues are a poor representation of 
s.p. energies; see [70, 71] for a recent study. In a DFT ap- 
proach, however, it is assumed that the generalized form 
of the energy density may effectively account for beyond 
mean-field effects such as particle- vibration couplings. In 
addition to this theoretical argument, let us recall that 
s.p. energies are not observables but model-dependent 
quantities extracted experimentally from binding ener- 
gies of excited states in odd nuclei. Systematic errors 
can thus be reduced by working exclusively with binding 
energies. 

To this end, we computed a number of onc-quasi- 
particle (q.p.) configurations for the odd-mass neighbors 
of 1^0, 40ca, 48Ca, s'^Ni, i32sn, and ^ospb. Calcula- 
tions were done at the equal filling approximation, which 
is an excellent approximation to the full time-reversal, 
symmetry- breaking blocking scheme [69]. Blocking q.p. 
states induces a small shape polarization [49], which in 
turn leads to a fragmentation of spherical s.p. orbitals of 
angular momentum j into 2j + 1 levels 51 = — j, . . . , + j. 
In principle, the "experimental" s.p. energy should be 
the average energy over all the 2j + 1 blocking configura- 
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FIG. 9: (Color online) Single-proton energies for ^"^Pb, cal- 
culated with unedfO and unedfI EDFs compared to the ex- 
perimental values of Ref. [72] . 



FIG. 8: (Color online) Two-neutron (top) and two-proton -48 
(bottom) separation energy residuals between unedfI results S^P ^^^^^ 
and experiment. 



tions. However, a state with the projection fl ean belong 
to any spherical orbit j > which could potentially 
complicate the identification for low-fi values. We there- 
fore associate a spherical orbit with spin j to the blocking 
configuration with the maximum projection ~ +j. We 
have verified that the amount of splitting in all these nu- 
clei does not exceed 150 keV: for the proton 113/2 orbit 
in ^°®Pb, it is 77 keV, and for the neutron I/7/2 state in 
'^^Ca, it is 124 keV. This implies that the error induced 
by cherry-picking a single state out of a 2j + 1 nonde- 
generate candidate instead of the average is only on the 
order of 60 keV. 

Figure 9 compares single-proton energies in ^ospj^ com- 
puted with UNEDfO and unedfI with the experimental 
values. Although the differences between these two func- 
tionals are small, we note that UNEDfI improves slightly 
the description of several high-j levels except for a too 
low position of the l/ig/2 proton orbit that reduces the 
size of the Z=82 proton gap. The situation is similar for 
the neutron single-particle energies in ^"^Pb. For both 
unedfO and unedfI the effective mass is close to one, 
which probably explains the fairly good reproduction of 
the level density in ^°*Pb. In the lighter doubly magic 
nuclei; differences between unedfG and unedfI s.p. en- 
ergies are somewhat larger. The magic gaps in Ca iso- 
topes are now better reproduced, although the iV = 28 



C. Superdeformed States and Fission Barriers 

Table IX lists the excitation energies of superdeformed 
(SD) fission isomers in the actinide region and SD band- 
heads in the mass A ~ 190 region calculated with the 
UNEDfO and unedfI parameterizations. Contrary to fis- 
sion isomers, SD bandheads in neutron-deficient lead and 
mercury isotopes were not included in the objective func- 
tion, since the prolate-oblate shape coexistence effects, 
not captured by current functionals, are well known in 
these nuclei [73, 74]. Indeed, calculations with unedfI 

predict an oblate ground state at /3=— 0.2 0.15 in all 

the Hg-Pb isotopes considered, coexisting with a slightly 
higher spherical minimum. By contrast, the ground state 
of the three lead isotopes is spherical with unedfO. The 
fact that the spherical configuration is disfavored in these 
nuclei can be traced back to a too-low Z~82 spherical 
proton gap in UNEDFl, (Sec. VB). 

All values listed in Table IX were obtained with the hf- 
BTHO code using the same large HO basis as used for the 
optimization. In particular, the deformation of the ba- 
sis was spherical for the ground state, and was deformed 
with /32 = 0.4 for the FI; see Sec. IIIB. As can be seen 
from Table IX, the optimization improves dramatically 
the rms deviation for the actinide nuclei included in the 
fit, going from 3.02 MeV in unedfO to 0.23 MeV in un- 
edfI. At the same time, the optimization deteriorates 
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TABLE IX: Excitation energies (in MeV) of fission isomers in 
the actinides and superdeformed bandheads in the neutron- 
deficient Hg and Pb nuclei calculated with hfbtho. The val- 
ues predicted with unedfO and unedfI are compared with 
experiment. 



Nucleus 


unedfO 


unedfI 


Exp. 


Ref. 


236u 


5.28 


2.42 


2.75 


[50] 


238 


5.73 


2.71 


2.557 


[50] 


240 


5.74 


2.51 


2.8 


[50] 


2«Cm 


5.27 


1.85 


1.9 


[50] 


^^^Hg 


6.33 


2.62 


5.3 


[75] 




7.27 


3.79 


6.017 


[76] 


192pj^ 


5.20 


1.25 


4.011 


[77] 


194p^^ 


5.99 


1.99 


4.643 


[78] 


196p^^ 


7.26 


3.52 


5.63 


[79] 



the description of SD excitations in the Hg and Pb iso- 
topes. To understand this behavior, we again compute 
the effective surface coefficient a^'^f (16) for ^^'^Pb. It 
is particularly low, 16.0 MeV, for unedfI. Indeed, it is 
17.64 MeV for unedfO, 17.1 MeV for SLy4, 16.50 MeV 
for BSk6, and 16.36 MeV for SkM*. In addition, the re- 
duced Z=82 magic gap in unedfI energetically favors 
deformed and SD states. Consequently, both bulk energy 
and shell effects of unedfI conspire to reduce the exci- 
tation energy of SD states in the Pb isotopes. In view of 
the major shape coexistence effects recalled earlier, this 
behavior is not too worrisome. 

In the A 190 region, the experimental uncertainty 
of the SD bandhead comes from the extrapolation of the 
rotational band down to spin 0"*". The associated error 
is estimated to be very small, around 5 keV. In the ac- 
tinides, experimental excitation energies of FI are usually 
determined with larger uncertainties. While the experi- 
mental error bar is only about 5-10 keV for 236,238u^ j^- 
grows to about 200 keV for ■^^^Pu and •^''^Cm. For the 
fission isomer of ^"'"Pu, recent measurements lower its 
excitation energy by about 500 keV to roughly 2.25 MeV 
[84] compared with the standard value [50] adopted in 
this work. Because of the relatively large uncertainty, 
Wi = 0.5 MeV, adopted for FI energies in our objective 
function, these experimental uncertainties are not going 
to significantly alter the final optimization. For future 
work, however, better experimental determination of FI 
bandheads should become a priority. 

Our long-term goal is to develop an EDF that can ac- 
curately predict and describe fission observables in heavy 
and superheavy nuclei. We present here some results of 
spontaneous fission pathway calculations with unedfI. 
All fission calculations were done with version 2.49t of the 
code HFODD [85] that can break all self-consistent sym- 
metries along the fission path. At each point along the 
collective trajectory, the HO basis deformation and fre- 
quencies are determined from a standard nuclear surface, 
parametrized by surface deformations a\^. The deforma- 



tions were chosen to minimize the energy as a function 
of the requested quadrupole moment Q20, according to 
the following expressions: 

3 2 

0^20 = A Q20 + B Q20 + C Q20, /I y\ 

a40 = 0.01, ^ ' 

with A = 3.16721 X lO-Sb-3, b = -2.75505 x lO'^b-^, 
C — 0.00954925 b~^, and all remaining deformations 
c^Xfj. = 0. The basis contains up to = 31 shells and 
up to A^ = 1140 states. Such an extended basis was 
previously applied in the systematic study of fission bar- 
riers in the transfermium region and yielded excellent 
convergence [8, 86]. We have checked that hfodd, with 
a reduced HO basis as compared with hfbtho calcula- 
tions, reproduces the hfbtho energies of FI bandheads 
in Table IX with an accuracy of 100-200 keV. We consider 
this agreement satisfactory considering other uncertain- 
ties involved. 

As an example. Fig. 10 displays the potential energy 
curve of ^^"^Pu as a function of the mass quadrupole mo- 
ment Q20- Triaxiality and reflection asymmetry effects 
are included for all calculations. The large-deformation 
behavior of the potential energy curve obtained by un- 
edfO is typical in this region, so that the outer barrier 
heights obtained by unedfO are systemmatically much 
higher than empirical values. We note that the un- 
edfI functional yields both the inner and outer barrier in 
^"^"Pu fairly close to experiment. Both unedfO and un- 
edfI functionals also yield the g.s. binding energy that 
is close to the empirical value. 




50 100 150 200 



Q2o(b) 

FIG. 10: (Color online) Fission pathway for ^''"Pu along the 
mass quadrupole moment Q20 calculated using hfodd with 
SkM*, unedfO, and unedfI EDFs. The ground state energies 
have been normalized to zero. En, Ea, and Eb indicate, 
respectively, the experimental energy of fission isomer and 
the inner and outer barrier heights [80]. 

The UNEDFl results for fission barrier heights in se- 
lected actinide nuclei are listed in Tables X (inner bar- 



13 



TABLE X: Empirical and theoretical inner barrier heights Ea (in MeV) for selected actinide nuclei. The rms deviations from 
experiment are AEa are shown in the last row. 



Nuclide 


Exp. [80] 


HFB+Fit [81] 


ETFSI [82] 


FRLDM [83] 


HFB-14 [31] 


SkM* 


UNEDFl 


236-[j 


5.00 


5.52 


5.20 


4.45 


5.52 


6.93 


6.39 


238-y 


6.30 


5.80 


5.70 


5.08 


5.93 


7.25 


6.50 


238py 


5.60 


5.57 


5.40 


5.27 


5.96 


7.39 


6.83 


240p^ 


6.05 


5.89 


5.80 


5.99 


6.49 


7.51 


6.77 


242p^ 


5.85 


6.02 


6.20 


6.42 


6.81 


7.44 


6.59 


244p^ 


5.70 




6.40 


6.59 


6.85 


7.82 


6.10 


2«Cm 


6.65 


6.20 


6.10 


6.56 


6.75 


8.76 


7.12 


244Cm 


6.18 


6.18 


6.40 


6.92 


7.10 


8.81 


6.99 


2«Cm 


6.00 


6.00 


6.50 


7.01 


7.31 


8.41 


6.69 


2«Cm 


5.80 




6.50 


6.80 


7.25 


7.94 


6.12 








0.47 


0.75 


0.87 


1.97 


0.79 



TABLE XI: Similar to Table X except for the outer barrier heights Eb (in MeV). 



Nuclide 


Exp. [80] 


HFB+Fit [81] 


ETFSI [82] 


FRLDM [83] 


HFB-14 [31] 


SkM* 


UNEDFl 


236-[j 


5.67 


6.03 


4.00 


5.03 


6.03 


6.70 


5.56 


238-y 


5.50 


6.17 


4.90 


5.64 


6.48 


7.36 


6.42 


238py 


5.10 


5.35 


2.90 


4.47 


5.24 


5.99 


4.62 


240p^ 


5.15 


5.73 


3.40 


4.91 


5.61 


6.40 


5.42 


242p^ 


5.05 


5.61 


3.60 


5.72 


6.02 


6.90 


6.20 


244p^ 


4.85 




3.90 


6.47 


6.25 


7.49 


6.50 


Cm 


5.00 


4.90 


1.70 


4.45 


4.51 


6.31 


4.08 


Cm 


5.10 


5.10 


2.10 


5.07 


4.83 


7.00 


5.03 


Cm 


4.80 


4.80 


2.40 


5.87 


5.23 


7.42 


5.51 


248 Cm 


4.80 




2.60 


6.65 


5.25 


7.32 


5.55 








2.11 


0.94 


0.70 


1.89 


0.84 



rier) and XI (outer barrier). For comparison, we also list 
the empirical barrier heights from the Reference Input 
Parameter Library (RIPL-3) [80] ; the HFB fission barri- 
ers obtained by fitting the neutron-induced fission cross 
section [81]; and predictions of ETFSI [82], FRLDM [83], 
and HFB-14 [31] models, together with hfodd calcula- 
tions with the SkM* EDF. (For SkM* predictions includ- 
ing the energy correction due to the rotational zero-point 
motion, see Ref. [6].) Overall, the description of experi- 
mental data by UNEDfI is very reasonable, with the rms 
deviations from experimental values of Ea and Eb com- 
parable to the values obtained in more phenomcnolog- 
ical models. One can thus conclude that fission barri- 
ers are reliably described at the HFB level with the un- 
edfI functional. This result is remarkable since it was 
obtained by adding only four excitation energies to the 
dataset. 

There seems to exist an interesting relation between 
barrier heights and the surface thickness. We have eval- 
uated the surface thickness in ^°*Pb from the charge form 
factor (Tch as defined, e.g., in Rcfs. [4, 23] and found 
= 0.932 fm for unedfO and (Tch = 0.907 fm for UN- 
edfI. This is to be compared with the measured value of 



fch = 0.913 fm [23]. It is apparent that the EDF which 
does well on fission barriers also performs well for surface 
thickness. The functionals SV-min and SV-bas which in- 
cluded (Tch in the fit, and yield values around 0.91 fm 
for this quantity, happen to perform well concerning bar- 
rier heights [48]. The relation between fission barriers 
and surface thickness deserves closer inspection in future 
work. 



D. Neutron Drops 

Recently, there has been a considerable interest in 
studies of inhomogencous neutron matter by consider- 
ing finite systems of neutrons; specifically, neutron 
drops [87, 88]. Since neutron drops arc not self-bound 
[89] , an external potential must be used to confine them. 
By studying neutron drops, one can test different ab ini- 
tio approaches and their correspondence to DFT calcu- 
lations [87, 88]; investigate the validity of the density 
matrix expansion [90] ; and develop a theoretical link be- 
tween neutron-rich nuclei and the neutron matter found 
in the neutron star crust [91]. 
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FIG. 11: (Color online) Comparison of unedfO and unedfI 
predictions for the energy A''-neutron drops trapped in a HO 
potential with tko = 5 MeV and 10 MeV with the AFDMC ab 
initio results of Ref . [881 . 



Figure 11 presents the results of unedfO and unedfI 
calculations for neutron drops confined by two external 
HO traps with fiw = 5 MeV and 10 MeV. The DFT re- 
sults are compared with ab initio AFDMC benchmark 
calculations of Ref. [88] employing the AV8' nucleon- 
nucleon and Urbana IX three-nucleon force. As can 
be seen, unedfO reproduces AFDMC results well, es- 
pecially since the functional has not been constrained to 
finite neutron matter. The agreement with unedfI cal- 
culations is also good, especially for a softer trap with 
huj = 5 MeV. 

In future EDF optimizations we shall include ab initio 
predictions for neutron drops into the dataset. By pro- 
viding unique constraints on finite neutron matter, such 
pseudo-data are expected to improve the description of 
very neutron-rich nuclei and diluted neutron matter. The 
results shown in Fig. 11 indicate that unedfO and un- 
edfI functionals represent excellent starting points for 
such optimizations. 



VI. CONCLUSIONS 

By performing nuclear energy density optimization at 
the deformed HFB level, we have arrived at the new 
Skyrme parameterization unedfI. Our main focus was 
to improve the description of fission properties of the 
actinide nuclei and to provide a high-quality functional 
for time-dependent applications involving heavy systems. 
The only notable change in the form of the energy den- 
sity as compared with our previous work [25] was the 
removal of the center-of-mass correction. For the x^" 
minimization, we used the derivative-free POUNDerS al- 
gorithm. Compared with unedfO, the dataset was en- 



larged by adding ground-state masses of three deformed 
actinide nuclei and excitation energies of fission isomers 
in 236,238u^ 240py^ ^j-^^j 242q^-^^ p^j. ^^iq optimal param- 
eter set, we carried out a sensitivity analysis to obtain 
information about the standard deviations and correla- 
tions among the parameters. We conclude that unedfI 
remains as robust under a change of individual data as 
unedfO. 

Overall, unedfI provides a description of global nu- 
clear properties that is almost as good as that of un- 
edfO. Not surprising, the quality of data reproduction is 
slightly degraded: by adding a new type of data (fission 
isomers), one is likely to worsen the fit for other observ- 
ablcs. The most striking feature of unedfI is its ability 
to reproduce the empirical fission barriers in the actinide 
region. We find it encouraging that, by including only a 
handful of fission isomer bandheads, deformation prop- 
erties of the functional seem well constrained. Another 
unanticipated property of unedfO and unedfI is their 
ability to reproduce ab initio results for trapped neutron 
drops. This is significant because such pseudo-data will 
be used in future EDF optimizations. 

In addition to imposing new constraints on neutron 
drops, in the next step we intend to improve the spec- 
troscopic quality of UNEDF functionals by considering 
the experimental data on spin-orbit splittings and shell 
gaps. We shall also improve the density dependence of 
the kinetic term by adding new constraints on giant res- 
onances. Meanwhile, the functional unedfI developed 
in this work will be the input of choice for microscopic 
studies of the nuclear fission process. 
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